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Changes in climate and atmospheric C0 2 and nitrogen (N) over the last several decades have induced 
significant effects on forest carbon (C) cycling. However, contributions of individual factors are largely 
unknown because of the lack of long observational data and the undifferentiating between intrinsic factors 
and external forces in current ecosystem models. Using over four decades (1956-2001) of forest inventory 
data at 3432 permanent samples in maritime and boreal regions of British Columbia (B.C.), Canada, growth 
enhancements were reconstructed and partitioned into contributions of climate, C0 2 and N after removal of 
age effects. We found that climate change contributed a particularly large amount (over 70%) of the 
accumulated growth enhancement, while the remaining was attributed to C0 2 and N, respectively. We 
suggest that climate warming is contributing a widespread growth enhancement in B.C.'s forests, but 
ecosystem models should consider C0 2 and N fertilization effects to fully explain inventory-based 
observations. 



Forest ecosystems that sequester carbon (C) from the atmosphere play an important role in Earth's C budget 
by offsetting the increase in atmospheric carbon dioxide (C0 2 ) caused by fossil fuel emissions and land-use 
change 1 . Net primary productivity (NPP) is strongly affected by many factors, including climate, atmo- 
spheric C0 2 and nitrogen (N) deposition, making it a challenge to reconstruct long-term forest NPP responses to 
these factors 2 " 6 . The separation of NPP into impacts of climate change, C0 2 and N is a prerequisite for future C 
accounting schemes that aim to establish a clear relationship between human activities and C uptake. 
Biogeochemical C cycle-climate models may be the best way to reconstruct long-term NPP as affected by 
environmental changes, projections of NPP responses to the changing future climate, and partition NPP varia- 
tions into contributions from various factors 7 . Though these models have uncertainties, they present the most 
feasible way to achieve these objectives as short-term C cycle simulation would not be reliable for forests 6 . To 
reduce large uncertainties in model simulations, observational data with long enough duration must be used to 
constrain and refine models. The problem, however, is the unavailability of these data to determine whether an 
observed response of NPP would be sustained through time. A more fundamental problem in most current 
analyses, which has not been given sufficient attention or even overlooked, is our inability to differentiate NPP 
responses to variations in intrinsic factors (e.g., stand development) and external forces (e.g., global change). This 
differentiation is important because stand development plays a significant role on forest NPP. For example, a 
commonly observed pattern of decline in NPP with stand age in mature forests may have implications on the 
decreased NPP response to elevated C0 2 8 . A critical step, therefore, is to remove the effect of stand development 
on NPP before we proceed to determine the contributions of climate change, C0 2 concentration and N deposition 
to NPP enhancement. 
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Figure 1 | Relationships between flux-tower derived estimates of annual 
net primary productivity (NPP) and InTEC- simulated NPP for three 
chronosequences Douglas-fir stands in B.C. DF00 (♦), DF88 (•) and 
DF49 ( ) represent flux sites of Douglas-fir stands harvested in 2000, 
1988, and 1949, respectively. Solid and dash lines are the regression and 
95% confidence levels for mean prediction, respectively. 



In this study, we first used an Integrated Terrestrial Ecosystem C- 
budget (InTEC) model which integrates stand age development with 
time to reproduce interannual variability of annual NPP at three 
chronosequences forest flux sites of British Columbia (B.C.), 
Canada. Then, InTEC model was recalibrated to reconstruct the 
long-term (1956-2001) stemwood growth enhancement (i.e., after 
removal effects of stand development from the absolute growth rates) 
revealed from repeated measurements at 3432 permanent inventory 
plots (PSPs) of forests in coastal (1966 plots) and interior (1466 plots) 
regions of B.C.. Contributions of climate, C0 2 fertilization and N 
deposition on accumulated growth enhancements were then quan- 
titatively determined by setting the factor in question constant in one 
simulation (e.g., C0 2 -based partition was estimated by detrending 
climate and nitrogen deposition used to drive the simulations). Our 
approach will provide useful guidance for future analyses of climate 
and atmospheric changes on global C sequestration. 

Results 

For three flux sites, the NPP values simulated using InTEC model 
agreed well with those derived from flux measurements, with coeffi- 
cients of determination (R 2 ) equal to 0.94 (p < 0.001) for the overall 
dataset and R 2 of 0.85 (p < 0.001), 0.81 (p = 0.002) and 0.62 (p = 
0.004) for DF00, DF88 and DF49, respectively (Fig. 1). Through this 
site level validation, we gained confidence that InTEC is a reliable 
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Figure 2 | Comparison of time series of InTEC simulations with observations from permanent inventory plots for maritime and boreal regions, (a) and 
(b): stemwood growth simulation of InTEC (■) and inventory data (•) (R 2 = 0.27, p < 0.001 and R 2 = 0.39, p < 0.001 for maritime and boreal 
plots, respectively), (c) and (d): residual growth of InTEC ( ■ ) and inventory data (•) (R 2 = 0.54, p < 0.001 and R 2 = 0.77, p < 0.001 for maritime and 
boreal plots, respectively). Error bars indicate (±) stand deviations. 
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Figure 3 | Temporal changes of baseline stemwood growth of (a) 
maritime (1966 plots) and (b) boreal (1466 plots) regions from 1956 to 
2001. 

model to link historical climate data with recent flux- tower observa- 
tions in reproducing interannual variability of NPP across forests of 
various stand ages in B.C. 

Age response function parameters for the InTEC simulations were 
recalibrated to match observations at PSPs (Supplementary 
Information). Inventory data suggested that stemwood growth 
(G_SW) of the coastal area showed a distinct pattern over the mea- 
surement period with a decrease during the middle 1960s and then 
recovered (—3.3 Mg C/ha/y) by the end of this decade (Fig. 2a). A 
gradual increase was observed since 1970 and G_SW reached 3.5 Mg 
C/ha/y at the beginning of 1990. In particular, a steep increase in 
G_SW was found for the next ten years with a value of 4.2 Mg C/ha/y 
for 1999, showing a rate of increase for this decade that was roughly 
seven times larger than that in the previous decades. The baseline 
stemwood growth (B_SW) generally exhibited a slow decrease over 
time, reflecting a small negative effect of stand aging on the stem- 
wood growth due to the increase in average forest age of the popu- 
lation of PSPs (Fig. 3). The residual growth (R_SW, i.e., after removal 
effect of stand development) showed the same pattern as G_SW and 
a tremendous growth enhancement was observed for the last decade 
(Fig. 2c). The interior area went through a quite similar growth 
pattern but with a lower G_SW and R_SW than that in the coastal 
B.C. area. 

Simulations from the InTEC model generally followed the pat- 
terns of inventory observations for both regions. For the coastal 
area, the pronounced decrease of G_SW in middle 1960s was well 
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captured. The main inconsistency for this region between the model 
and measurement occurred in the last decade when a rapid increase 
of G_SW was observed in the inventory data while InTEC- simulated 
increase was less rapid. However, if the effect of stand age was 
removed, this inconsistency can be mostly reduced, as illustrated 
in Fig. 4c with comparable slopes of growth enhancement curves 
between InTEC simulation and inventory data. For the interior 
region, model results were also consistent with the inventory data 
in terms of capturing the slight increase of G_SW in middle 1970s, 
the quick increase after middle 1990s, and the general pattern of 
R_SW for the whole measurement period. 

Based on these evaluations of InTEC with two different datasets, 
we then partitioned growth enhancement into contributions from 
different factors including climate, C0 2 and N, and further separated 
the total climate effect into those from temperature, precipitation, 
radiation and water vapor pressure (Fig. 4). We found that climate 
was the dominating factor controlling on accumulated R_SW with 
70.9% and 79.4% for coastal and interior regions, respectively with 
the remaining signal attributed to C0 2 and N. In the coastal region, 
18.7% of accumulated R_SW was attributed to C0 2 and 12.8% for N. 
Similar results were observed for the interior region that C0 2 and N 
contributed 15.8% and 6.9% to R_SW, respectively. Partitioning the 
climate effect into components due to temperature, precipitation, 
radiation and water vapor pressure suggested that temperature is 
the main factor (48.4% and 57.6% for coastal and interior regions, 
respectively). Contributions of other three factors differ by some 
extent between these two regions. For the coastal area, contributions 
for precipitation, radiation and water vapor pressure were 20.2%, 
15.8% and 14.6%, respectively, while these values were 12.7%, 5.9% 
and 22.5% for the interior region. 

Discussion 

The most distinct feature of the InTEC model compared with other 
coupled C-climate models is its incorporation of the NPP-age rela- 
tionship of forests, by which InTEC considers the structural changes 
associated with stand development through time. This feature gives 
InTEC a capacity to remove the influence of the age factor on stem- 
wood growth, making it possible for the consequent analyses of the 
contributions from individual non-disturbance forcings. The 
removal of the age effect in this process is the necessary prerequisite 
since increasing evidence shows that stand age plays an important 
role in forest growth and neglecting the age effect may lead to incor- 
rect responses of NPP or tree growth to other factors 910 . Results of 
our analysis also supported this theory that by removal of baseline 
growth resulted from stand ages, correlations of residual growth 
between simulations and observations were improved substantially 
for both regions. 

Using over three thousand permanent sample plots over four 
decades, we observed a growth enhancement of forests in British 
Columbia, Canada. Contributions from climate change, C0 2 and 
N deposition were factored out with a process-based model. 
Climate change has been demonstrated as the main contributor to 
the growth enhancement in both coastal and interior regions, 
explaining over 70% of the trend. This dominant role of climate 
change may indicate that global warming has markedly affected these 
high-latitude ecosystems. Such a pronounced impact of climate 
change on C dynamics is not unusual 2,411 . Our analysis further 
showed that the temperature is the most important factor of climate 
change in the high-latitudinal regions in B.C., and the underlying 
mechanism is the significantly increased growing season length 
(GSL) for two regions over the observation periods (Fig. S2). 

C0 2 and N generally contributed similar fractions (below 20%) to 
the overall growth enhancement (see Fig. S3 for data of C0 2 and N) 
in the two regions studied here, yet they are also important to com- 
pletely explain the growth enhancement. Direct comparison of the 
C0 2 contribution with previous analyses might not be appropriate 
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Figure 4 | Partition of the growth enhancement observed from inventory into contributions from climate (Climate), C0 2 fertilization (C0 2 ) and 
nitrogen deposition (N) for the (a) maritime and (b) boreal regions. Partition of growth enhancement due to climate into contributions from 
temperature (tern), precipitation (pre), radiation (rad) and water vapor pressure (vap) for the (a) maritime and (b) boreal regions. 



since the C0 2 fertilization effect is reported to be not uniform 
through the various forest ecosystems and could vary with stand 
characteristics 7,812 " 14 . While there is no unique sensitivity of NPP 
to increasing C0 2 concentration, the identification of its contri- 
bution is still important for regional C dynamics. The impact of N 
deposition would depend on whether the ecosystem is N limited and 
for how long it sustains. In N limited forest ecosystems, the likely 
effect of N deposition is the increased foliar N concentration with a 
positive effect on photosynthetic rates, which is observed for oceanic 
spruce stands 11 . Similar results were also observed by other studies 
showing enhanced production of leaves and wood in response to N 
deposition 15 . A simulation indicates that the N deposition generally 
increased NPP by 17% of US forests during last century 12 . Our simu- 
lations generally showed slightly lower (12.8% and 6.9% for coastal 
and interior regions, respectively) N contributions. However, these 
data still indicates that N contribution should be considered when 
explaining the growth enhancement in this region. 

The contributions of both C0 2 and N to growth enhancement 
were slightly higher in the coastal region than in the interior region. 
Mechanistic analyses of the differences in these contributions might 
be difficult at present, but a possible reason could be the interactions 
of these factors with water availability. The reason is that the effect of 
C0 2 fertilization is balanced between NPP stimulated by increased 
C0 2 and NPP decreased by potential drought. Therefore, even at 
elevated C0 2 , NPP may still decrease under severe drought so that 
the C0 2 fertilization would be more evident in humid environments 
compared to arid or semi-arid regions 16 . The coastal region has more 
annual precipitation (1942 mm ±174 mm) than the interior region 



(673 mm ± 64 mm) so that forests in this region endure little 
drought and can realize the full potential of the C0 2 fertilization 
effect, while forests in the interior BC may experience periodic sum- 
mer droughts that limit the duration of C0 2 fertilization 17 . It should 
be noticed that the sum of contributions from all factors are not 
100%, given that results from the portioning approach are approxi- 
mately because of interactions of various factors and the nonlinear 
model descriptions. However, the interaction was low (within 2.5%), 
probably due to low values of these variables (high values tend to 
reach the saturation in sensitivity) and relatively short time duration 
of observation. 

The NPP -age curve is the key component of the InTEC model and 
is consequently vital for the main conclusions of our analysis. Here 
we provide an additional discussion on this relationship to support 
our results of partition. Stand age is an important regulator of annual 
NPP, which increases at the initial several years of growth and max- 
imizes at the mature age. After that, annual NPP generally shows a 
decrease trend and finally stabilizes in a certain level to maintain C 
neutral for old forests. The Douglas-fir stands of the Pacific 
Northwest are of high productivity (flux sites in this study) and 
annual NPP increases to around 1200 g C/m 2 for DF49 which has 
an age of 65 (Figure 1). However, the PSPs used had an average age 
above 70 (Figure SI) and therefore, evidences should be provided to 
support the decreasing NPP response to an age above 70. For this 
purpose, we substantially reviewed the NPP-age response of forests. 
Previous analysis showed that annual NPP generally maximizes 
around an age of 35-60 for three chronosequences forests using 
ground measured data 18 . A later analysis also indicated that annual 
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NPP would decrease after 40 years of growth for most North 
American forests 19 . Recently, a study further illustrated that the max- 
imum annual NPP would be observed at an age around 55 for 
Canadian forests 20 , which agrees with results that shows the NPP 
of Douglas-fir stands of the Pacific Northwest generally enter into 
the decreasing trends above 70 21 . The currently most comprehensive 
analysis on NPP-age relationship of North American forests demon- 
strates that the timing of maximum annual NPP would fall in the 
range between 20-60 using forest inventory data covering more than 
150,000 permanent field sample locations 10 . With these previous 
analyses on NPP-age relationships, in particularly those in Canada 
and Pacific Northwest regions, we suggest that our conclusions of 
partition based on NPP-age curve of InTEC would be reliable and 
acceptable. 

With flux measurements and over four decades of inventory data, 
we partitioned stemwood growth enhancement of forest in B.C. into 
contributions of climate (temperature, precipitation, radiation, 
humidity) and atmospheric chemistry (C0 2 concentration, nitrogen 
deposition) using the process-based InTEC model that considers 
stand development in simulating processes. Our results suggest that 
climate warming is contributing a widespread and significant growth 
enhancement in B.C.'s forests, but that inclusion of C0 2 and N 
fertilization effects in the model are also necessary to fully explain 
the observed growth enhancement. Among the various climate fac- 
tors, temperature was found to be the main factor inducing the 
growth enhancement in these high-latitude regions. 

Methods 

Flux data and permanent inventory plots. Two types of data were used in this study, 
including continuous flux measurements at three chronosequences forest sites in B.C. 
and 3432 permanent inventory plots (PSPs) of forests in coastal (1966 plots) and 
interior (1466 plots) regions of B.C. (Figure SI). Plots for these two regions belong to 
maritime and boreal ecozones, respectively. Detailed descriptions of the flux sites 
(Table SI) and PSPs were provided in Supplementary Information. 

Data processing. Half-hourly ecosystem C0 2 flux data were continuously measured 
at each site with the eddy-covariance technique 22 and these data were obtained from 
Fluxnet- Canada archives (http://www.fluxnet-canada.ca). Detailed procedures for 
flux data processing (gap-filling and partition net ecosystem exchange (NEE) into 
gross primary productivity (GPP) and total ecosystem respiration (Re)) were 
provided in the Supplementary Information. 

Repeated measurements of tree diameter and height between 1956 and 2001 were 
provided by the Inventory Branch of the BC Ministry of Forests, Lands and Natural 
Resource Operations and used to calculate stand-level stemwood biomass using 
allometric equations 23 ' 24 . Mean annual stemwood biomass growth (G_SW) was then 
calculated over the five- or ten-year measurement intervals. The combined effect of 
intrinsic factors on tree growth, herein referred to as baseline stemwood growth 
(B_SW), was modeled by fitting Weibull distribution functions to the age response of 
observations for each type of dominant species 25 . The remaining (unexplained) 
variation of growth (R_SW) was then calculated by subtracting B_SW from G_SW. 
More specific descriptions on these inventory-based variables are provided in the 
Supplementary Information. 

Integrated terrestrial ecosystem C-budget (InTEC) model and its inputs. The 

InTEC model is a process-based carbon cycle model 18 . Different datasets were used in 
the simulation of NPP at flux towers and PSPs (See Supplementary Information). To 
bridge NPP simulations and stemwood growth observations, we assumed that 
stemwood growth consists 3 1 percent of NPP for coniferous stands in InTEC 
simulations. The partition of stemwood growth enhancement into contributions 
from climate, C0 2 and N factors was realized by setting the factor in question 
constant in one simulation 26 . Detailed introduction of the InTEC model and its inputs 
are provided in Supplementary Information. 
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